The 4 He tetramer ground state in the Faddeev-Yakubovsky 
differential equations formalism 



I.N. Filikhin, S.L. Yakovlev 
Dept. of Math. & Comput. Physics, Institute for Physics 
Sankt - Petersburg State University 198504 Ul'yanovskaya 1, Petrodvorets, St. Petersburg, Russia 

V.A. Roudnev 

Physics Department, University of South Africa, P.O. Box 392, Pretoria 0003 

B. Vlahovic 

Physics Department, North Carolina Central University 
1801 Fayetteville St. Durham 27707 NC, and 
Jefferson Lab, 1200 Jefferson Avenue, Newport News, VA 23606, USA 

(February 2, 2008) 

Abstract 

The characteristics of the four 4 He atom cluster are investigated using the 
differential equations for Yakubovsky components. Binding energy, mean- 
square radius and density function are calculated for the ground state. The 
spatial properties of the cluster and its subsystems are studied. 
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I. INTRODUCTION 



In the present work clusters of three and four helium atoms are investigated. In recent 
years such systems had attracted growing attention due to development of experimental 
techniques allowing to observe and to study two and three-particle Helium clustersBi and as 
a result of theoretical developments which made possible the construction of realistic pair 
potentials ab initiJfin. 

The last two decades have been significant progress in the investigation of few-body 
systems. This progress rests on the development of new techniques for solving Faddeev 
and Yakubovsky (FY) equationsi'01ll0 for wave function components and growing power 
of computer facilities that allowed to reach higher accuracy in the branch of Monte- Carlo 
techniques. Though equivalent to Schrodinger equation the system of FY equations has 
significant advantages over the Schrodinger one. These advantages have their origin in 
a proper choice of wave function decomposition into the components. The equations are 
constructed so that only one pair potential enters the equation for a particular component 
provided that the interaction in the system is given in terms of pairwise potentials. This 
nice feature leads to very simple boundary conditions for the components and simplifies 
their numerical approximation significantly. Traditionally Yakubovsky equations are used 
in nuclear physicsEOHJ. In spite of strong mathematical background and effective numerical 
techniques developed for Faddeev and Yakubovsky equations only a few papers devoted to 
molecular physics exploiting the technique of Faddeev or Yakubovsky equations are known. 
Most of them are devoted to a system of Helium trimerEliHil. No four-body calculations of 
molecular systems based on Yakubovsky equations are known to the authors up to now. 
The aim of the present paper is to take the first step towards wider exploiting the rigorous 
and effective technique of Yakubovsky equations in molecular and chemical physics. 

The numerical scheme exploited in this paper to solve differential Yakubovsky equations 
is the cluster reduction method (CRM)lll which opens a way to construct small subspaces 
containing solutions of the equations. This method developed and repeatedly applied to the 
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problems of nuclear physics has been employed in the present paper to perform a calculation 
of bound state characteristics of a system of four helium atoms. The method allows to obtain 
not only the estimations for the binding energy of a four particle system but also the wave 
function of the system using computer resources economically. Another important advantage 
of CRM is applicability of the method to model-free calculations of multichannel reactions in 
a system of four particleJMMl. Being applied to systems of four helium atoms it can be used 
to calculate the reaction rates for He2 + He2 — > He3 + He and He2 + He2 — > He^ + He processes. 
Calculating s-wave model of He4 tetramer in this paper we make a step towards investigation 
of reactions in systems of four atoms on the base of strongly grounded mathematically correct 
methods. 

This paper contains three additional sections and a conclusion. In section 2 the model 
is described and the main equations are given. In the section 3 the method of solution is 
presented. Section 4 contains the results of numerical calculations. 



II. FORMALISM 

The calculations presented here were performed in the framework of differential 
Yakubovsky equations (DYE) for four particles in configuration space. The formalism of dif- 
ferential Yakubovsky equations was developed by S.P.Merkuriev and S.L.Yakovlev§. Here 
we only give a brief description of these equations touching upon the approximations that 
we use and emphasizing the advantages of using the DYE instead of the Shrodinger equation 
in investigations of bound states of four particles. Detailed and sound description of the 
equations and asymptotic boundary conditions can be found in the original paper0 and in 
more detail in the book by Merkuriev and Faddeevi. 



When considering a system of four particles it is convenient to use Jacobi coordinates. 
For identical particles there are two types of Jacobi coordinates, which correspond to differ- 
ent partitions of the four particle system into subsystems. The first type corresponds to the 
partitioning into a three-particle cluster and one separate particle (3+1 type). The second 
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one corresponds to the partitioning into two-particle clusters (2+2 type). The explicit ex- 
pressions for the Jacobi coordinates through the particle coordinate vectors r k , k — 1, 2, 3, 4 
are given by the formulas 

x 2 ,i = r 2 - ri , 

y3,2i = (ri + r 2 )/2 - r 3 , (1) 
24,321 = (ri + r 2 + r 3 )/3 - r 4 , 
for the first type of coordinates and 

x l, 2 = r 2 " r l i 

y 3 , 4 = r 4 - r 3 , (2) 

z 34 ,i 2 = (ri + r 2 )/2 - (r 3 + r 4 )/2 , 
for the second one. The Jacobi coordinates for partitions of the same type but with different 
distribution of particles among the clusters can be obtained by cyclic permutations of the 
subscripts enumerating the particles in (|l|) and (@). 

Suppose the Hamiltonian of a system of four particles has the form 

H = H + Y / V(x a ) , (3) 

a 

where H is the Hamiltonian of the system of free particles, a stands for two-particle sub- 
systems of the four-body system and V(x a ) is the potential of the interaction in the pair 
with index a. For the systems of identical particles the wave function \l/ can be expressed 
in terms of two Yakubovsky components U k , k—1,2. U 1 corresponds to the partition 3+1 
and U 2 corresponds to the partition 2+2. The expression of the total wave function of the 
system in terms of Yakubovsky components reads@lll 

¥ = (/ + P+ + P+P+ + P~)(I + P 4 + + p 4 -)c/ x + 

(4) 

+(/ + + Pf )(I + P + P + )U 2 . 

Here P + (P~) are the operators of cyclic (anticyclic) permutations of four particles, Pj 
correspond to cyclic permutations of three particles with fixed z-th particle. The Yakubovsky 
components U 1 and U 2 satisfy the following set of the equations^: 



{H + V{x) — E) U l + V(x)(P+ + P 4 ) U 1 = 

= - v(x) Up? + p + ) u l + (p+ + p+) p 2 



(5) 

(P + V{y) -E)U 2 + V(x)(P+P+) P 2 = 

= - y(x)(P + + P+)P + P 1 . 

Here we have omitted the subsystem index a in the notation of coordinates x since all the 
particles are identical. The advantages of using the equations for Yakubovsky components 
instead of solving directly the Shrodinger equation come from better localization of the in- 
teraction in configuration space. In the special case of identical particles the DYE can be 
written in term of the interaction potential of only one pair. As a result the numerical 
approximation of the Yakubovsky components is a much easier problem than the approxi- 
mation of the wave function. Detailed discussion of Yakubovsky equations can be found in 
the monograph by Merkuriev and Faddeevi. 

We solve the equations for Yakubovsky components in the s-wave approximation in which 
the angular momenta of the system of four atoms and all its subsystems are set to zero. The 
s-wave equations for Yakubovsky components U k , k—1,2 have the following form! 

1 

/ xy 
dv U (xx, yi, Z\) = 

1 1 1 

= —-v(x) du dv ( U 1 (x 2 ,y2,Z2) + 

' x 2 y 2 Z2 

+ ^^U\x^z,)), (6) 
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(hi + v(x) — e)U 2 (x,y,z) + v (x) U 2 (y, x, z) = 
= --v(x) / du U (x 4 , 2/4,^4) 

2 J XaVaZa 
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where x = |x|, y = \y\, z 



^ = -(^ + ^ 2 + ^ 2 ), 

hl = -(dl + dl + \dll 
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v(x) is s-wave component of the pair potential V(x). The coordinates Xj, yi, Zi, i =1,2,3,4 in 
the kernels of the equations (||) are defined by the following relations: 

X!=\—+y+xyv\ , y x = I (-x) + — - -xyv I , 



x 2 = xi, x 3 = xi, x A = y, 




III. METHOD OF SOLUTION 

The differential equations for the Yakubovsky components (§) are solved using the clus- 
ter reduction method (CRM). This method has been developed and applied beforellMzlll! 
to calculate the characteristics of bound states and low-energy scattering of systems of 
three and four particles. The cluster reduction method reduces considerably the computa- 
tional difficulties when solving DYE numerically. In the framework of the CRM Faddeev 
(Yakubovsky) components are decomposed in terms of the eigenfunctions of the Hamiltoni- 
ans of two (three) particles subsystems. As a result of the projection onto the elements of 
a biorthogonal basis we obtain the set of equations corresponding to the relative motion of 
clusters. A brief summary of the CRM from Yakovlev and FilikhinEl is given below. The 
Yakubovsky components U l , i—1,2 are written in the following form 

oo 

U i (x,y,z)=Y / <P}(^y)Fi(z), i = l,2. (7) 

1=0 

In the Eq. (^) the basic functions are the solutions of s-wave Faddeev equations for 
subsystems of types 3 + 1 (i = 1) and 2 + 2 (i = 2): 
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{~d 2 x - \d 2 y + v(x)}<ft(x, y) + v(x) [ dv^-tfix!, Vl ) = 

{-d 2 -d 2 y + v(x)}(f) 2 (x, y) + v(x)<ft(y, x) = e l 2 (j) 2 (x, y) . 

The set of functions {ipj} biorthogonal to the set {<$]} consists of the eigenfunctions of 
the equations adjoint to the Eq. (H) 

{-d 2 x -p 2 y + v(x)}^l(x,y)+ fdv-^vixMixuyx) = 

= e[^}{x,y) , ( 9 ) 
{-d 2 x -d 2 + v(x)}ip?(x, y) + v(y)ipf(y, x) = e l 2 ip?(x, y). 

A biorthogonal basis is required because the Faddeev operator is not self-adjointEl. Sub- 
stituting (|7|) into the Eq. @ and projecting onto conjugated basis {ipj} we obtain the set 
of integro-differential equations for the functions F{(z), describing the relative motion of 
clusters 

{-ld 2 z +e[-e}F l \z) = 

, oo / 11, 

= -5 E Ui(x,yMx) J du J dv{^ti(x 2 ,y 2 )F£(z 2 ) + 

k=0 \ -1 -1 y 

+ if£#(*s,vs)J?(*)}), ( 10 ) 

= - fc E y)\v(x) J dv^tKxi, y,)Fl{ Zi )\ . 

In these equations (.|.) means the integration over the variables x and y. The functions 
.F7(z) must vanish when z — > oo 

2?(*)~0, i = l,2, / = l,2,...,oo. 

The number of equations in the set depends on the number of the terms retained in the 
expansion of the Yakubovsky components, Eq. (|7]). Due to completeness of the set of the 
basic functions only a finite number N of such terms needs to be taken into account to 
support a stable numerical solution. 
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IV. RESULTS 



The solution of the Eq. (10) has been computed in the region f2 of the configuration 



space defined by the parameters R x , R y , R z : 

Q = {x, y,z : x < R x , y < R y , z < R z }. 

The values of these parameters were chosen to be R x =R y =R z =50 A. All the calculations were 
performed using model potentials HFDHE2I and HFD-b!. From the one hand according; 
to contemporary point of view these potentials give lower and upper limits for two-body 
binding energies correspondingly, and from another hand some four-body results for these 
potential models are known in literature^. 

The basic functions (pl(x,y), ipl(x,y), i = 1,2, I = 1,2, ...,N were calculated using 
the CRJVllll. To confirm the accuracy of the basic functions we were checking the basis for 
orthogonality using the condition 

((0j|^)-M<lO- 3 , i = l,2. 

The function <p\(x, y) for I = 1 and k = 1 in the Eq. (|7]) is the s-wave Faddeev component 
of the ground state wave function of the He3 system (trimer). The binding energy of the 
He3 ground state when computed for the HFDHE2 and the HFD-B potentials has the values 
-0.105 K and -0.118 K. Compared to the values reported by Carbonell et al.il and Kolganova 
et al.il our trimer is slightly overbound. 

The He4 (tetramer) binding energy, which has been computed using the potentials 
HFDHE2 and HFD-B, is given in the Tab. 1. In the same table we also quote the re- 
sult reported by Nakaichi-Maeda and LimEI. These authors used the formalism of integral 
AGS equationsil. In addition we include the results of the calculations of the mean square 
radius (< r 2 > l l 2 ) of the system, the mean square distance between Helium atoms, and the 
probability of forming the subsystem with cluster structure He3 + He. This last probability 
has been computed as 
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P 3 +i =<^Mi|* >, (11) 

where \I/ is the total wave function of the system, ipi is the ground state wave function of 
He3 and F\ is the function describing the motion of He3 trimer relative to a single He atom. 
As can be seen from the table, the contribution of the He3 + He state to the total wave 
function is considerable. 

The fast convergence of the cluster decomposition, Eq. fl7|) indicates the existence of 
clusters in subsystems. Particularly, in the 3+1 subsystem one needs to take into account 
two terms to get a stable binding energy. The binding energy computed by taking into 
account only the 3+1 component has the value -0.25 K, which is in good agreement with 
the value -0.24 K by Nakaichi-Maeda and Lim@. To achieve a binding energy calculation 
that is stable at the scale 10~ 2 K one needs to take into account six terms in the Eq. (^) 
which use components of both types (3+1) and (2+2). 

The analysis of the results of our computation of the characteristics of He3 and He4 
systems enables us to draw an analogy with the nuclear cluster systems 3a and 4a. Here 
the symbol a denotes a 4 He nucleus. The bound states of these systems correspond to 
the ground states of nuclei 12 C and 16 0, respectively. It was known (for exampldlll) that 
in these cases the three-body systems have no well defined clusters in their subsystems. 
However, in the four-body systems the cluster of 3+1 type is dominant i.e. in this case it 
is possible to separate a closely bounded cluster of three particles and fourth particle. The 
mean square radius of the nuclear systems 3a (< r 2 > x / 2 =2.33 fm) and 4a (< r 2 > 1 / 2 =2.54 
fm)0 increases with the number of particles. The situation is similar for the system under 
consideration. In particular for the potential HFDHE2 (HFDH-B) the mean square radius 
of trimer is 6.7 A(6.5 A) and that of the tetramer is 7.4 A(6.9 A) for the same potentials. 

To characterize the spatial distribution of the particles constituting the tetramer we 
have computed its wave function (|j) for the HFDHE2 potential. The density function p(r) 

oo 

is depicted in Fig. 1. normalized with the usual condition / p(r)dr = 1. 

o 

To study the spatial position of the Helium atoms in the tetramer we have used the total 
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wave function of the system. For comparison analogous calculations were performed for the 
trimer. The He 4 (He 3 ) wave function depends on six (three) variables. These are moduli 
of the Jacobi coordinates x, y, z (x, y) and cosines of the angles between vectors x, y, z 
(x and y) u — v = w = The most probable configurations of the relative 

position of the particles forming the He 4 (He 3 ) system was calculated as the coordinates of 
the maximum of the square of the total wave function. For He 3 system we found x=3.6 A, 
y—S.l A, u—0, and for He 4 system x=3.6 A, y—3.1 A, z=2.9 A, u—0, v—0, w—0. These 
configurations are shown in Fig. 2. For the ground state of the trimer the center of the 
He atoms masses arrange themselves at the vertices of the equilateral triangle with sides 
as large as 3.6 A(Fig. 2. a). For the ground state of the tetramer the three Helium atoms 
are located at the vertices of an equilateral triangle with sides as large as 3.6 A, while the 
most probable position of the fourth Helium atom is at a distance of 2.9 A in the direction 
perpendicular to the plane of the three particle system (Fig. 2. b) and through the center of 
the equilateral triangle formed by them. One should not be surprised by the predominance of 
the tetrahedron configuration if one takes into consideration the identity of particles and the 
s-wave approach that has been used for the description of the tetramer. One can compare the 
location of density function maximums with the positions of potential energy minimums. 
These positions differ noticeably that demonstrate the essentially quantum nature of the 
system. The minimums of potential forms a configuration of equilateral tetrahedron with 
the side of 3.0 Awhereas the maximums of density function are located on the vertices of 
tetrahedron with the side of 3.6 A. 



V. CONCLUSION 

By applying the method of cluster reduction we have solved numerically the s-wave 
differential equations for the Yakubovsky components for a system with four 4 He atoms. 
Binding energy, mean-square radius and density function are calculated for the ground 
state. The results of the calculations are in good agreement with those of Nakaichi-Maeda 
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and which were performed using the integral equations. The configurations with He3 

cluster and separated helium atom dominates in the He4 cluster. This behavior is analogous 
to that of the nuclear 4a particles system0. The most probable spatial configuration of the 
four Helium atom system is the tetrahedron with sides as large as 3.6 A. 
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FIGURES 

FIG. 1. The comparison of the probability densities for three- and four He atom system: solid 
curve corresponds to the ground state of the He4 system, dashed curve corresponds to the ground 
state of the He3 system (HFDHE2 potential). 

FIG. 2. The most probable configurations of Helium atoms: a) the ground state of the He3 
system, b) the ground state of the He4 system. The figures show the numbers of atoms and 
distances between their centers of mass (HFDHE2 potential). 



14 



TABLES 



TABLE I. 4 He4 tetramer binding energy (£4), mean-square radius (< r 2 > 1 ^ 2 ), mean square 
distance between Helium atoms (< x 2 > 1 ^ 2 ) and the contribution of cluster subsystems of the He3 
+ He form (P3 + i). 



Potential 




HFDHE2 


HFD-B 




present work 


S. Nakaichi-Maeda and T.K. Lim@ 


present work 


E 4 , K 


-0.39 


-0.394 


-0.41 


< r 2 > 1 /2, A 


7.4 




6.9 


< X 2 > 1 /2, A 


11.1 




10.3 


P3+I 


0.75 




0.81 
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